Thermodynamics of water modeled using ab initio simulations 



Valery WebeiQ 

Physical Chemistry Institute, University of Zurich, 8057 Zurich, Switzerland 

D. Asthagirj^ 

Department of Chemical and Biomolecular Engineering, 
Johns Hopkins University, Baltimore, MD 21218 
(Dated: August 9, 2010) 

We regularize the potential distribution framework to calculate the excess free energy of liquid 
water simulated with the BLYP-D density functional. The calculated free energy is in fair agreement 
with experiments but the excess internal energy and hence also the excess entropy are not. Our work 
emphasizes the importance of thermodynamic characterization in assessing the quality of electron 
density functionals in describing liquid water and hydration phenomena. 



The liquid-vapor coexistence for water is thoroughly 
characterized experimentally, and the excess free energy 
(chemical potential) of a water molecule in the liquid rel- 
ative to the vapor, /i^, a basic descriptor of the liquid- 
vapor equilibrium, is very well-established. The chemi- 
cal potential, together with its derivatives, especially the 
temperature derivative (the excess entropy) , are essential 
quantities in understanding hydration phenomena and 
chemical transformations in the liquid state. 

For ah initio simulations of water, two earlier studies 
[Tl[2] have sought /x^. The first study revealed how the 
overbinding of water by PW91 (PBE) functionals leads to 
an overly negative chemical potential and ultimately an 
over-structured fluid. The second explored the coexisting 
densities of the liquid and vapor phases P] and found a 
higher than normal vapor density, indicating overbinding 
of molecules by the BLYP functional. Their limitations 
notwithstanding, these early studies were incisive in eval- 
uating the fluid simulated by ah initio dynamics. 

In this Letter we present a rigorous calculation of the 
excess free energy using a theoretical approach that has 
matured over the last few years [SHS], one that renders 
calculating free energies from ab initio simulations far 
less daunting than before. Together with an independent 
calculation of the excess internal energy, we obtain the 
excess entropy as well. A key finding of this work is that a 
satisfactory agreement with experiments of the chemical 
potential can mask errors in the excess energy and the 
entropy, quantities that are crucial in understanding the 
thermodynamics of hydration. 

The relation between /i™ and intcrmolecular interac- 
tions is given by the potential distribution theorem [7] 

I3fil^ = \nj eP'Pie)de, (1) 

where e — Un ~ Un-i — CAv is the binding energy of 
the distinguished water molecule with the rest of the 
fluid. Un is the potential energy of the iV-particle sys- 
tem at a particular configuration, Un~i is the potential 
energy of the configuration but with the distinguished 
water removed, and is the potential energy of the 



distinguished water molecule solely. P{e) is the probabil- 
ity density distribution of e and is obtained by sampling 
many configurations of the system, is the excess free 
energy in the liquid relative to an ideal gas at the same 
density and temperature. Based on the experimental co- 
existence densities [8] at 298 K, ^™ = —6.3 kcal/mol. 

A naive application of Eq. [T] to liquid water will fail 
because the high energy regions of P{e), reflecting the 
short-range repulsive interactions [H IHl [TO], are never 
well sampled in a simulation. We resolve this difficulty 
by regularizing Eq. [l] [10,. Consider a hard-core solute 
of radius A centered on a distinguished water molecule: 
the hard-core solute excludes the remaining water oxygen 
atoms from within the sphere of radius A. The chem- 
ical potential of the hard-core solute, ^ffg, is assumed 
known. With this construction, Eq. [l] can be rewritten 
as [HEIIIO]: 

/J/i^'^ = /3Mhc + Ina^o + In j P{e\n). - Q)e^'de . (2) 

Here P{e\n\ — 0) is the binding energy distribution of the 
distinguished water molecule conditioned on there being 
no {n\ = 0) water oxygen atoms within A of the dis- 
tinguished oxygen atom. By moving the boundary away 
from the distinguished water, we temper the interaction 
of the distinguished water with the rest of the fluid; in- 
deed, for a sufficiently large A, P{e\n\ = 0) is expected to 
be well-described by a Gaussian [4] . The fraction of con- 
figurations that do not have any oxygen atoms within A 
of the distinguished oxygen is Xq; this is the n = mem- 
ber of the set {x„} of coordination states sampled by the 
distinguished water molecule and characterizes the inter- 
actions between the distinguished water and the contents 
of the inner-shell. The excess chemical potential of the 
hard-core solute Pfi^Q = — Inpo: where po is the fraction 
of configurations in which a cavity of radius A is found 
in the liquid, pq is the n = member of the set {pn} 
of occupation numbers of a cavity and characterizes the 
packing interactions in hydration llj. 

Since a;o and po are directly obtained from the simula- 
tion record, we only need explicit energy evaluations to 
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compute the outer-shell contribution 



(3) 



Here we calculate the outer-shell contribution using an 
alternative expression 



In P'-°\e\nx = 0)e-^'de, (4) 



where P^^'^{e\n\ = 0) is the binding energy distribu- 
tion in the uncoupled ensemble [6l [THj; that is, after 
we find a suitable cavity in the liquid, we insert a test- 
particle in that cavity and assess P'^^^e\nx = 0). The 
superscript (0) indicates that the test-particle and the 
fluid are thermally uncoupled, ft is advantageous to 
use Eq. |4] over Eq. |3] because the uncoupled distribution 
P^^\s\n\ — 0) has a higher entropy than the coupled 
distribution P{e\nx — 0), rendering the calculation of 
the free energy of inserting the particle (Eq. |4| more ro- 
bust than extracting it (Eq. [3| [12J. Further, when xq is 
small, it is difficult to characterize P{e\n\ = 0) reliabily, 
a limitation that does not apply to P^'^^(e|nA = 0). 

For A sufficiently large such that P*^°-'(e|nA = 0) is 
Gaussian distributed with a mean {e\nx = 0)o and vari- 
ance (fc^jriA = 0)o (the subscript emphasizes that 
the test-particle and the fluid are thermally uncoupled), 
Eq. [4] becomes 



/3M™tcr = (el^A = 0)o - 
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(fc2|nA=0)o. (5) 



We apply the above framework to water simulated with 
the BLYP-D electron density functional. (In BLYP-D, 
the BLYP [ISllM] functional is augmented with empirical 
corrections for dispersion interactions [H].) We simulate 
the liquid at a density of 0.997 g/cm^ (number density of 
33.33 nm~'^) and temperature of 350 K. The higher tem- 
perature effectively weakens the bonding [S] and has been 
found necessary to model the real liquid at the standard 
density and T = 298 K [TB] . 

We simulate water with the Cp2k code^Tj and in the 
NVE and NVT ensembles; the electronic structure cal- 
culations are exactly as in our earlier study [5], and only 
salient differences are noted here. For the studies at 
NVT, we use the hybrid Monte Carlo method and ex- 
tended the simulations in our earlier study of a 32 water 
system [S]. In the NVE ensemble, we simulated systems 
with 32 and 64 water molecules. 

The production phase of the hybrid Monte Carlo lasted 
2585 sweeps («130 ps). The system temperature was 
350 K. Simulations in the NVE ensemble lasted 200 ps 
of which the last 170 ps were used for analysis. The 
initial configuration for the 32 particle system was ob- 
tained from the end-point of our earlier study [S]. The 
initial configuration for the 64 water system was obtained 
from an equilibrated configuration of SPC/E water |18J 



molecules. In the production phase, the average temper- 
ature was 357 ± 26 K (362 ± 19 K) for the 32 (64) water 
system. The energy drift was less than 1.8 K (1.2 K) for 
the 32 (64) water system. 

To calculate Pouter' construct a Cartesian grid 
within the simulation cell. The oxygen population within 
a radius A of each grid center is calculated. In instances 
where there are no {n\ = 0) water molecules, we insert 
a test water molecule in the cavity, assess its binding 
energy, and thus construct P''^\e\n\ = 0). (We draw 
water molecules from the simulation cell, randomly ro- 
tate it about an axis through the oxygen center, and use 
the resulting configuration as the test particle.) 

In the NVE simulations, we first find approximately 
3500 cavities from configurations sampled every 50 fs. 
(The target count of cavities is met by adjusting the 
spacing of the Cartesian grid between approximately 0.3 
to 0.96 A.) Thus over 110k(« 32 * 3500) binding energy 
calculations are used to construct P'^°^(e|nA = 0). In 
the NVT simulations, the grid spacing for finding cavi- 
ties to compute P^'^\£\n\ = 0) was fixed at 2.0 A. For 
A = 3.0 A this gave 98 cavities. In this instance, five dif- 
ferent orientations were used for each test particle giving 
about 16k(w 5 * 32 * 98) binding energies to construct 
P''^\£\n\ = 0). Subsequently, we resampled the trajec- 
tory using a finer grid to better estimate the probability 
of finding po- Throughout, uncertainties in xq, po, and 
{e~^'^\n\ = 0) were estimated using the Friedberg and 
Cameron [19] block transformation procedure [20] . 

Figure 0shows P'-°'){e\nx ~ 0) for different cavity radii 
for the 32 water system. As the figures shows, the low 
energy region is well-characterized for inner-shell radii 
considered here; and it is this low-e part of the distribu- 
tion that is most important for calculating /ioutcr (Eq- 14| . 

Figure [T] shows that the agreement between 
P^^Helnx = 0) obtained in the NVE and NVT 
simulations is satisfactory, although there is somewhat 
more scatter in the low-energy region for results with 
NVT. This is a consequence of using an order of 
magnitude less data (16k versus 110k) to construct 
p(o)(e|nA -0). 

The probability xq is a direct measure of the solute- 
water interactions within the inner-shell, and perhaps not 
surprisingly, for large A, this quantity is very low. In 
this case, one can model {x„} with a maximum entropy 
approach [51] [52] using the robust estimates of mean and 
variance of {xn} and thus secure xq. But for large A 
[31 [S] , the underlying assumption of Gaussian occupancy 
statistics may not be valid. Here, we instead use Bayes's 
theorem in the form [3] 

p(l|n + 1) 



— Pn-\-l 



Em>lPmP(l|w) 



(6) 



where p(l|m) is the conditional probability of finding one 
water molecule at the center of the cavity given that m 
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FIG. 1. Distribution of interaction energies P''''(e|nA = 0) 
for a distinguished water molecule in the center of an empty 
cavity of radius A. The system comprises 32 water molecules. 
The red triangles and blue circles are from the NVE simu- 
lation. The solid line defines the Gaussian model for the re- 
spective distribution. The black dots are from the NVT sim- 
ulation. The A = 3.3 A curve has been shifted downwards by 
3 units for clarity. Observe that the thermodynamically im- 
portant (Eq. |4| low-energy tail is well-characterized. Further, 
for A = 3.3 A, both wings of the distribution are tempered 
and the Gaussian model becomes a good approximation. 

water molecules are present in the cavity. For this pur- 
pose, following an earlier study the center is any point 
within 0.15 A of the geometric center of the cavity. 

Figure [2] depicts {x„} for the 64 particles system. The 
fit using Eq. |6] is in excellent agreement with the actual 
data. This gives us confidence in the estimated xq for 
A — 3.3 A (Fig. [2]), an instance where xq was not observed 
in the simulation. As Fig.[2]shows, the maximum entropy 
model can lead to errors in k^ThiXQ by about a kcal/mol 
or more, especially for the large A. 

Table |l] collects the various components of the excess 
chemical potential. Across the range of A, the calculated 
values are within a kcal/mol (often less) of each other and 
the experimental value of —6.3 kcal/mol for liquid water 
at 298 K. Comparing the 32 and 64 water results shows 
modest system size effect of about k-oT. This difference 
primarily arises from a less favorable /^outor for larger 
system: for a large cavity in a small simulation cell, the 
outer-shell water molecules pack tightly at the interface 
and lead to the lower /^outer- 

We can calculate the excess entropy of hydration per 
particle from excess chemical potential and the mean 
binding energy, (e), of a distinguished water molecule 
in the fully coupled simulation. First consider the ex- 
cess internal energy. For the 32 (64) water system 
(e) = -28 kcal/mol (-27.3 kcal/mol). Thus the heat 
of vaporization per particle AiJvap ~ — (^)/2 + k^T — 
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FIG. 2. Observed and estimated {a;„} for a A'^ = 64 water 
molecule system and different inner-shell radii A. Symbols 
denote the observed data. The solid line is obtained using 
Bayes' theorem (Eq. |6|. The dashed lines are maximum en- 
tropy fits using the mean and variance of {xn} and a flat prior 
[22j . Observe that {a;„} obtained using Eq. [6] well-describes 
the simulation results. 

14.6 kcal/mol (14.2 kcal/mol) is substantially in error 
relative to the experimental value of about 10.5 kcal/mol 
at 300 K 0. 

Neglecting small contributions due to the thermal ex- 
pansion and compressibility of the medium, the excess 
entropy per particle d [10] is S'^'^/nfes « {e)/2kBT- 
fi^^/kBT. Based on the estimates of (e) and above, 
it is clear that TS°^ for BLYP-D water is also in error. 

Several important lessons emerge from the present 
study. First, while BLYP-D at 350 K appears to ade- 
quately describe the structure [5] and the hydration free 
energy relative to experiments, the agreement comes due 
to balancing errors in the excess entropy and excess in- 
ternal energy. Second, while dispersion interactions are 
regarded as necessary in describing water, the empirical 
dispersion correction in BLYP-D overbinds the material 
and the temperature of 350 K does not adequately com- 
pensate for this overbinding. Finally, the BLYP-D func- 
tional at 350 K should not be expected to describe the 
hydration thermodynamics of aqueous solutes, especially 
properties such as the excess entropy of hydration, and 
thus also transport properties of solutes in the medium. 
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TABLE I. Chemical (kBTlnxo), packing (fcerinpo), and long-range interaction {fil^ter) contributions to the excess chemical 
potential (/i™) of water obtained from NVE molecular dynamics simulations and NVT hybrid Monte Carlo simulations 0. 
The radius of the inner shell is A (in A). The average temperatures of the = 32 and = 64 water systems are, respectively, 
357 ± 26 K and 362 ± 19 K. The temperature is 350 K for simulation in the NVT ensemble. Where available, the inner-shell 
chemical contribution is evaluated directly from the data, else the estimate using E g. [6] is used. The outer-shell contribution 
is calculated by numerical integration of Eq. |4] estimates from a Gaussian model [5| (data in parenthesis) are provided for 
comparison. Energies are in kcal/mol. Uncertainties are at la level. 
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